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SUMMARY 


A quadrilateral shell element, CQUAD4*, has been added to level 15.5 
and subsequently to level 16.0 of NASTRAN. The element exhibits doubly- 
curved surfaces and uses bi-quadratic interpolation functions. Reduced 
integration techniques are used to improve the performance of the element 
in thin-shell problems. Several details of previous authors' (ref. 1) 
work are clarified with respect to the present NASTRAN implementation. 

The creation of several new bulk data items is discussed along with a 
special module, GPNORM, to process SHLNORM bulk data cards. In addition 
to the theoretical basis for the element stiffness matrix, consistent mass 
and load matrices are presented. 

Several potential sources of degenerate behavior of the element are 
investigated. Guidelines for proper use of the element are suggested. 
Performance of the element on several widely-published classical examples 
is demonstrated. The results show a significant improvement over pre- 
sently available NASTRAN shell elements for even the coarsest meshes. 
Potential applications to two classes of practical problems are discussed. 


INTRODUCTION 


Until recently, only the CQUAD2 and its analog CTRIA2 were available 
in NASTRAN for analyzing shells of arbitrary geometry. Compared to 
current shell element technology, these elements are subject to the 
following limitations: 

o Faceted (flat) surface geometry is poorly adapted to model 
curved shapes. 


* After the initial implementation of the new element was completed, the 
authors became aware of a similar proprietary element under develop- 
ment by the MacNeal-Schwendler Corporation which used the name CQUAD4. 
The reader should take care not to confuse these two identically 
named elements, since it is our understanding that the formulation 
and performance are quite different. 
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o Lower order polynomials used in membrane formulation cause 
element to be excessively stiff for in-plane deformation. 

o While the enforced linear variation of the normal slope along 
the sides of the element guarantees interelement compatibility, 
it causes the bending behavior of the element to be quite 
stiff as well. 

o In problems exhibiting thick shell and/or three-dimensional 
behavior over certain regions, the CQUAD2 element is an 
inadequate model and is difficult to interface with three 
dimensional elements. 

To alleviate these problems development work on the present (CQUAD4) 
element was begun with the intention of implementing it in NASTRAN 
Level 15.5. The efforts were partially successful but full implementation 
was not achieved until NASTRAN Level 16.0 became available last year. 

The choice of the element was primarily influenced by the need to 
accurately represent curved surfaces as well as thick shel 1/3-D behavior. 
Such extremely accurate elements as Cowper's (ref. 2) and Dupuis' 

(ref. 3) were rejected due to the present authors' preference to adhere 
to the standard six degrees of freedom (dof) preferred by the majority 
of the user community. Although the theoretical development has often 
been presented elsewhere (refs. 1, 4, 5, 6, 7, and 8), we choose to 
repeat enough of the development to clarify certain issues which caused 
difficulty in the present implementation. 
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SYMBOLS (Scalars) 


Values are given in both SI and U.S. Customary Units. The 
measurements and calculations were made in U.S. Customary 
Units. 


a, b Plate edge dimensions 

0 Shell flexural rigidity [Et 3 /12(l~v2)) 


E 


Elastic modulus 


g li’ g 2i’ 
g 3i ’ h 3i 
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N 1 (I, n. O 

P 

P 

q 

R 

t. 

i 

U, V, w 


U 1 , V 1 , w 1 


u i 1 

v t ■ w t 


i 

i 

i 

v lx> 

iy’ 

v lz 

i 

i 

i 

v 2x> 

V 2y ’ 

v 2z 

v 1 

v 3x J 

v 1 * 

3y> 

V 1 

v 3z 


Components of interpolation derivative arrays 

Shear correction factor 
Structure length dimension 
Interpolation function for node i 
Constant pressure load on element 
Concentrated load magnitude 
Displacement value 
Mean (midsurface) radius 
Thickness at node i 

Translational displacements at a point in basic coordinate 
system 

Translational displacements at a point in local coordinate 
system 

Translational displacements at node i 

Components of unit vector defining local x axis at node i 
Components of unit vector defining local y axis at node i 
Components of the shell normal at node i 
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x, y, z 


x' , y' , z' 

V V z i 

a i* p i 


Y^ly! . Y X t Z l 1 

Y I I 

V Z 


e x'x' * e y'y' ’ 

£ Z'2' 


Basic cartesian coordinate variables 
Local cartesian coordinate variables 
Basic cartesian coordinates at node i 
Rotational displacements at node i 

Shearing strain components in local coordinate system 
Direct strain components in local coordinate system 


v 

u) 

i, n. C 
4-f , Dj » ^ 

CT x'x' ’ a y'y' * 
°z' z 1 

x x 1 y * * x x' z' ’ 

x y' z' 

(b . . 

V 1 J 

P 


Poisson's ratio 

Natural frequency of structure 

Curvilinear coordinate variables 

Curvilinear coordinate at node i 

Direct stress components in local coordinate system 

Shearing stress components in local coordinate system 

Component of transformation matrix 
Mass density per unit volume 
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SYMBOLS (Matrices and Vectors) 


[ALPHAT] 

[BETAT] 

[B 1 ] 

IB-*] 

[DELTAT] 

[D 1 ] 

f 

CG-1 

[H-] 

[J] 

[K] 

[M] 
n 

[N] 

s 
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-> 
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s- 





V 


li J 


-> 


2i 




3i 


Diagonal matrix of nodal thicknesses times local x rotation 

Diagonal matrix of nodal thicknesses times local y rotation 

Strain-displacement relation referenced to local coordinates 

Strain-displacement relation pertaining to node i 

Array of translational displacements at nodes 

Constitutive relation in local coordinate system 

Consistent load vector for element 

Derivative array transformed to local coordinates 

Derivative array transformed to local cocrdi nates per- 
taining to 3/84 operator 

Jacobian matrix relating (x, y, z) and (4, g, 4) systems 

Element stiffness matrix referenced to basic coordinates 

Element mass matrix referenced to basic coordinates 

Third row of Jacobian - Interpolated value of nodal normals 

Vector of nodal interpolation functions 

Array of nodal interpolation functions 

First row of Jacobian - vector tangent to surface 
C = const 

Second row of Jacobian - vector tangent to surface 
4 = const 

Unit vector tangent to surface 4 - const, defining local x 1 
axi s 

Unit vector tangent to surface 4 = const, defining local y 1 
axis 

Unit vector normal to surface 4 = const, defining local z 1 
axis 

Unit vectors defining local tangent coordinates at node i 
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[V1TAN] 

[V2TAN] 

[V3N0RM1 

[XCOORD] 

h 

t 

P 

[ 6 ] 

S' 

w 

[n] 


Vectors defining the coordinate system for nodal rotations 

Array of local x direction vectors at nodes 

Array of local y direction vectors at nodes 

Array of shell normals at nodes 

Ar.’iy of nodal coordinates 

Displacement dof at node i 

Collection of nodal displacement vectors 

Strain components in local coordinates 

Local/global transformation matrix (direction cosines) 

Stress components in local coordinates 

Transformation from (x 1 , y' , z') system to (4, n> O 
system 

Differential operator matrix for computing strains 
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THE STIFFNESS MATRIX 


Basic Assumptions 

Figure 1 shows the geometry of a typical element. The curvilinear 
coordinate system (£, n » 0 is used where £ and n lie in the middle 
surface of the element while £ is directed through the thickness. Each 
of these coordinates is allowed to vary from -1 to +1 on opposite faces 
of the element. We adopt the customary assumption of shell theory that 
the strain component (e 2 , z i) in the thickness direction is neglected 

compared to the other strains. The input items describing the element 
geometry include the basic coordinates at each of the eight mid-surface 
nodes (GRID cards) plus the vectors normal to that surface at each node 
(SHLNORM cards). The length of each normal vector is taken to be the 
thickness at that node. The thickness is interpolated quad“atically over 
the element. At present only homogeneous, isotropic, materials (MAT1 
cards) are allowed. The element is not available for heat transfer 
problems nor are thermal load vectors calculated. 


Interpolation Functions 

The nodal coordinates are related to the basic coordinates by the 
equation: 
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( 1 ) 


Details of the biquadratic interpolation functions (Ni) and their deriva- 
tives are given in Appendix A. 
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The translational displacements are chosen as u, v, and w in the x, 
y, and z directions respectively, Two rotations, ot. and p., are defined 

^ ^ * 

about the local axes, v^.. and v^. , tangent to the mid-surface at each 
node (i). The choice of these two local axes is discussed in Appendix B. 

We may relate the nodal displacements to the continuous displace- 
ment representation in a manner analagous to equation (1). 


where 
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[BETAT] 


hh 0 


o p 2 t 2 . . 


. . o 

. . 0 


o .... p n t f 


Coordinate Transformations 

Tfve relation between the curvilinear coordinates (|, q, £) and the 
basic coordinates (x, y, z) is commonly called the Jacobian, [J] , defined 
as: 


9x/9£ 

9y/94 

9z/9£ 

9x/9q 

9y/9q 

9z/9q 

9x/9£ 

3y/9t 

9z/9£ 


Using equation (X) we can write out these expressions as: 

{ 9x/94) 

s = = [XCOORD] • + | ■ [V3N0RM] • afi/8£ (4) 

az/H 


i ax/aq 

ay/9n 


= [XCOORD] • 9t5/9q + | • [V3N0Kf‘(„ • 9t$/9q (5) 


where 


9z/9q 


j 3x/9C ) 
\ 3y/sC j 

( 9*/3£ ) 
9^/9| = 


\ • [V3N0RM] • H 


< aN-j/94 3N 2 /94 


afi/Oq = < 3N 1 /9n 9N 2 /9q 


. 9Ng/94 > 

. 9Ng/9n > 


The explicit forms of 9Ni/9£ and 9Ni/9q are given in Appendix A. 

Physically the s and t vectors may be considered tangent to the surface 
4 = constant, while the n vector is merely the interpolated value of the 
node normals and may not be exactly normal to that surface at the position 

(4, n, C). 
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Perhaps the most confusing point of the cited references is the use 
of still another local coordinate system for definition of the stresses 
and strains. The need for this additional (x 1 , y l , z‘) system arises 
from the definition of the basic shell assumptions (particularly the 
neglect of the through- the- thickness direct strain, c z 1 2 , ) . We wish to 

define, at any point in the element, a local (z 1 ) axis which is normal 
to the surface £ = constant along with two other orthogonal axes (x 1 , 
y') which are tangent to that surface. Since we have previously deter- 
mined that s and t are tangents to the surface, we can determine a 
normal vector as: 



= s x t 

(7) 

and 

V a $ /|$ 1 

n n 1 n 1 

(8) 


The other two unit vectors defining the local axes (v and v^) are 

computed in a manner analagous to that given in Appendix B. Thus, x 1 is 

measured along the v_ vector, y' is measured along the v. vector, and z' 

* ^ 

is measured along the v vector. We can define the transformation [0] 

as: n 


So 



(9) 

( 10 ) 


and 



(ID 


The Strain-Displacement Relation 


Using the newly developed local system and noting that e , , is 


neglected, we write the basic definition: 
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( 12 ) 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


We make use of equations (10) and (3) to develop the relation 

’ a/8x 1 ' 

3 /ay 1 
,a/3z\ 


(a/3x) (a/an 

= [ef 1 l a/ay > = [a] -1 [J]" 1 < a/an > 
(a/az) \b/K) 


(13) 


Define 


[<l>] = [0] T Uf 1 


(14) 


Where [0] 1 is equal to [0]^, since [0] is defined to be an orthonormal 
matrix. Then [()>] will have the form: 


11 

*12 

21 

*22 

31 

*32 


33 


(15) 


A complete derivation of the terms in [(J>) is given in Appendix C. 

However, [J] ^ is best evaluated numerically at each integration point 
and cannot be written out explicitly. We combine equations (13), (14), 
and (15) with (12) to arrive at: 
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(16) 


Finally we use (11) along with (2) to substitute the appropriate expres- 
sion for the displacements < u 1 v' w 1 > 


s' = [fi] [0] T [DELTAT] • t + [fl| • [0] 
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By carrying out the indicated operations to allow ^he differential 
operator [fi] to appropriately interact with £ and N and by rearranging 
terms, we arrive at the relation: 

e' = [B']t (18) 

where ^ ^2 ... Sg > ^ 

*c» 

and t. = < u . v. w. a. B. > 

1 T 1 1 1 r 1 

The explicit form of [B 1 ] is shown in Appendix D. 

The Stress-Strain Re 1 a t i o n 


Again referring to the primed local coordinates, the constitutive 


law is: 

P = [0*] e 1 

(19) 

where 
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and (for 

a homogeneous isotropic material): 
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here k is used to improve the shear representation. The displacement 
assumption causes the shear to be constant through the thickness, whereas 
the proper distribution is closer to parabolic. The ratio of the strain 
energies of the two distributions (parabolic/constant) is 1.2 which is 
substituted for k. 

The Element Stiffness Matrix (Subroutine KQUAD4) 

The standard virtual work arguments lead to the stiffness computation 
as follows: 

[K] = J [B‘] T [D 1 ] [B 1 ] dVol (21) 

Vol 

The usual volumetric measure is dx dy dz. Here the variables of integra- 
tion are £, r\ t and £. The conversion of the cartesian volume to the 
curvilinear volume is via the Jacobian. Thus, 
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dVol = dx dy dz - det [J] d£ dr| d£ 


( 22 ) 


111 

So [K] = J f f [B'] T [D 1 ] [B 1 ] det [J] d£ d n d£ (23) 

-1 -1 -1 

The integration is carried out numerically using two Gauss points in 
each direction. While capable of properly integrating the element's 
volume, this "reduced" integration is not sufficient to exactly evaluate 
the complex polynomials produced by equation (23). This implies that, 
while ultimate convergence is assured, the behavior will not be either 
bounded or monotonic. However, several authors (ref, 4 and 9) have 
shovjn that, by purposely underestimating the energy, the performance of 
the element is enhanced. By taking [J] and [0] to be invariant through 
the thickness, it is possible to explicitly carry out the integration in 
£. We choose not to do so, however, in order to ensure complete gen- 
erality of the formulation for both thin and thick shell cases. 

Stress Recovery (Subroutines SQUD41 and SQUD42) 

Once the elements are assembled and the system equations solved for 
the displacements, the user needs to know the element stresses as well. 
Combining equations (18) and (19) with o now known, we obtain; 

a' = [D 1 ] [B ' ] 3 (24) 

Recall however that, in general, [B 1 ] is a function of the curvilinear 

coordinates (£, Hi £)• <?' is therefore also a function of these coor- 

dinates, so we must choose which points we will use for stress evalua- 
tion. It is known that the numerical integration points are the best 
"samples" of the overall element stress field. Unfortunately the values 
of £ = ±0.57735 do not give the maximum stresses through the thickness 
if bending is present. We have compromised to select the eight points 
given by £ = ±0.57735; q = +0.57735, and £ = ±1.0 to allow evaluation of 
the stresses at the top and bottom surfaces of the element (c.f. diagram 
in Appendix A). Since the values of cf^ « x i > a y'y' anc * x x'y'* are eva ^ ua ' tec * 

in the (x 1 , y 1 , z') local system, the stress directions may not be mean- 
ingful to the user. Consequently, the principal stresses (ct^, o^, and 

t 1 are also calculated and form the additional portion of each line 

tTiaX 

of output. 


THE CONSISTENT MASS MATRIX (Subroutine MCQUD4) 


For simplicity we will neglect the rotational inertias associated 
with the a and p degrees of freedom. This assumption is particularly 

^Notice that and x , are zero on the top and bottom surfaces. 
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appropriate for thin and moderately thick shells. It allows us to 
reference everything to the mid-surface (£ = 0). In particular: 


[M] = X p [N] T [N] dVol (25) 

Vol 

Choosing p - constant and making use of the above assumption (i.e. 

1 

X d£ S t) 

-1 

[Ml a p x t [N] T [N] dA (26) 
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In general p may be allowed to vary quadratically over the element in a 
manner similar to the thickness. This feature is not required for most 
cases of interest. 


Since £ = 0 on the midsurface, we must resort to the device of 
computing the unit area in curvilinear coordinates as (using equations 
(4) and (5)): 

dA = dx dy = |s x t| d| dq (26) 

where js x t| may be interpreted as the projection on the normal vector 
v^ of the normal vector associated with infinitesimal area, dx * dy. 

Thus, 

11 _ -I- — 

[M] = p X X t [N] 1 [N] |v | di dq (27) 

-1 -1 n 

In this case the full three-point Gauss integration must be used to 
properly evaluate the expression. 


THE CONSISTENT LOAD VECTOR (Subroutines PL0AD4 and PWORK) 


We derive the expression for a constant pressure (p) normal to the 
mid-surface of the element. As before, a quadratic variation of the 
pressure would cause no inherent difficulties. The development is 
entirely analagous to that used for the consistent mass matrix. Thus, 
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Again, the three-point Gauss rule is used to evaluate the expression. 


SPECIAL NASTRAN CONSIDERATIONS 


New Bulk Data Cards 


Three new Bulk Data cards have been added to NASTRAN in conjunction 
with the new element. They are: 

CQUAD4 “ describing the element connectivity 
PL0AD4 - specifying the elements to which constant pressure (p) 
is applied at the mid-surface 

SHLNORM - inputting the direction vector of the normal to the shell 
surface at each grid point. 

A complete description of each of these items is found in Appendix E. 
Module to Process 5110:" Normals * 

A new module, GPNORM, has been coded which converts the “external" 
grid point ID's on a SHLNORM card to the appropriate internal SIL's. 

The module also transforms the normal vector into the basic coordinate 
system for the problem and writes the results on the output data block 
SHLNRM. The DMAP calling sequence for the module is: 

GPNORM G E0M1 , EQ EXI N , BGPDT , CSTM / SHLNRM $ 

GPNORM must be added to the DMAP rigid format immediately prior to the 
TA1 module. SHLNRM must be added as the final input data block of TA1. 

Augmented ECPT and EST Data Blocks '*' 

The make up of the EST (and by analogy the ECPT) for the CQUAD4 
element follows the standard format for the first 43 words. 

Word Contents 

1 Element ID 

2 Material ID 

3-10 8 Grid Point SIL's 

11-42 8 sets of Grid Point CSID's 

plus basic x, y, z coordinates 
43 Element Temperature 


*The idea to use GPNORM to process the shell normals as well as the 
technique for augmenting the ECPT and EST data blocks is credited to 
Miles Hurwitz of NSRDC. 



The formulation of the element requires the components of the shell 
normals. These must be appended to the EST by subroutine TA1A. 

44-67 8 sets of x, y, z components of 

the shell normals in basic 
coordinates 

Subroutine TA1B performs a similar augmenting process for the ECPT data 
block. 


VERIFICATION 


Consider the limiting case of a square, simply supported flat plate* 
subjected to two load conditions, 1) a central load normal to the plate 
and 2) a uniform normal pressure. Figure 2 indicates that excellent 
convergence to the Timoshenko (ref. 10) result can be obtained with a 
1 x 1 or at most a 2 x 2 grid. Note, however, that the usual bound 
theorems are not available with this particular element due to the use 
of reduced integration. 

The next step in analytical complexity is . -presented, by the prob- 
lem portrayed in Figure 3, a pinched cylinder with free ends,* Accord- 
ing to Timoshenko, the radial deflection at the point of application of 
the load, for the geometry given, should be -2.76 mm (-0.1087 in.). 
Timoshenko's result is based on an assumption of inextensional deforma- 
tion which neglects the middle surface strain of the shell. The CQUAD4 
element gives a slightly higher result of -2.89 mm (-0.1139 in.) for a 
325 degree of freedom model of one-eighth of the cylinder. Cantin and 
Clough (ref. 11) predict a deflection of -2.87 mm (-0.1128 in.) using a 
cylindrical shell element model with 1200 degrees of freedom for one- 
eighth of the cylinder. Therefore, the CQUAD4 element, although not 
monotonic in convergence, does give excellent results for a minimum 
number of degrees of freedom. 

An example problem which has become a classic for checking the 
response of shell -type elements is shown in Figure 4. The example is a 
cylindrical shell roof loaded by its own weight.* The ends of the shell 
are supported by diaphragms and the sides are free. It should be noted 
that two "exact" solutions have been quoted by various researchers. 

These two solutions may be attributed to Scordelis and Lo (ref. 12) and 
Cowper, Lindberg, and Olson (ref. 2). 

Scordelis and Lo based their calculations on the theory of Gibson 
(ref. 13) essentially using shallow shell equations. Cowper, Lindberg, 
and Olson claimed that the shallow shell approximations were not used 
consistently when particular loadings were considered. They expanded 
the trigonometric representation of the load variation up to second 
order within each element by means of a Taylor Series. In addition 

*When proper symmetry conditions are applied, only 1/4 or 1/8 of the 
entire structure need be modeled in each of these cases. 
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Cowper, et al. performed the integration of both the stiffness and load 
matrices over the area of the actual shell surface. Hence, the primary 
difference between the two "exact" solutions is the manner in which the 
consistent load matrices are formulated. Ihe formulation of the present 
element follows more closely the tethud of Scordelis and to and hence 
will be compared to their "exact" solute 

Table 1 gives the computed displacements based on the grids defined 
in Figure 4. Note that excellent convergence is obtained. Figure 5 
compares the predictions of the CQUAD4 element with the CQUAD2 element 
as well as slightly different formulations of the CQUAD4 element in the 
MARC (ref. 14) and SUPERB (ref. 15) finite element programs. Based on 
the results indicated in Figure 5 the CQUAD4 element is judged to be the 
most accurate. 

To demonstrate the applicability of the CQUAD4 element in modeling 
dynamics problems, consider the rectangular cantilever plate vibration 
problem reported by Zienkiewicz (ref. 16) (see Figure 6). Compared in 
Figure 6 are test results by Plunkett (ref. 16) and finite element pre- 
dictions based on a non-conforming triangle by Zienkiewicz and results 
from the CQUAD4 element. Note that even the two element idealization 
with the CQUAD4 element gives excellent results for the first four 
modes. 


POTENTIAL SOURCES OF ELEMENT DEGENERACY 


Three potential sources of element degeneracy were investigated. 

The first, non-rectangul arity of the mesh, is illustrated in Figure 7. 

One quarter of a simply supported flat plate subjected to uniform pres- 
sure was modeled as shown with angular "offsets" or variations in the 
mesh rectangul arity of up to 30°. As indicated in Figure 7 by the 
displacement prediction for the center of the plate, variations of up to 
20° resulted in only 2 % variation in deflection compared to the regular 
rectangular grid. The 30° variation resulted in a 7% difference. It 
would appear from these results that, for most applications, non- 
rectangularity will not have a significant effect on the results. How- 
ever, care should be taken to maintain small angle variations of less 
than 30° as good practice. 

The second potential source of degeneracy investigated was the shell 
thinness ratio, t/R. A pinched cylinder example was once again selected 
and the t/R ratio varied from 0.1 to 0.0001 as shown in Figure 8. By 
examining the product of the radial deflection at the point of load 
application and the flexural rigidity of the shell it is evident that no 
numerical instability exists even for very thin shells. Notice that for 
thicker shells (e.g. t/R > 0.1) Timoshenko's assumption of thin shell 
behavior is increasingly violated and some deviation of the finite 
element results from the classical solution is obtained. 
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The third potential source of element degeneracy to be investigated 
was similar to the mesh non-rectangularity of the flat plate problem. 
This example considers the event of a misalignment of the edges of the 
element with the directions of principal curvature in a shell idealiza- 
tion. Such a discretization may necessarily occur as a result of com- 
plex intersections of several shell elements. The pinched cylinder 
problem discussed previously was chosen as a simple limiting case, A 
2x2 grid was selected for one-eighth of the cylinder and symmetry 
conditions were enforced. The element edges in the circumferential 
direction were allowed to vary from the direction of curvature as shown 
in Figure 9. The authors found that the radial deflection at the point 
of load application was virtually unaffected for the cases examined. 

It may be concluded that, based upon the previously described 
investigations regarding element degeneracy which could possibly result 
from potential misuse, the CQUAD4 element appears to be exceptionally 
stable. Care should be used, however, in maintaining "relatively" 
rectangular element configurations. 


APPLICATIONS 


At least two potential sources of application for the CQUAD4 ele- 
ment exist in the offshore industry. Offshore drilling and production 
platforms are typically either a space frame of tubular members, commonly 
called a steel jacket structure, or a reinforced, prestressed concrete 
structure, commonly called a gravity structure. The welded intersec- 
tions of tubular members in a steel jacket are called tubular joints and 
represent sources of potential fatigue problems due to high stress 
concentrations. The CQUAD4 element represents a significant increase in 
computational accuracy compared to the CQUAD2 element for conducting 
stress analyses of these tubular intersections. 

The reasons for the improved accuracy are two-fold. The curved 
surface of the CQUAD4 element is its most obvious advantage. It was 
often necessary to use excessive CQUAD2 elements in otherwise coarse 
mesh regions of the joint model just to approximate the cylindrical 
geometry. An extremely important but less obvious advantage of the 
CQUAD4 element is its higher-order representation of the displacements, 
strains, and stresses, without having to expend any additional degrees 
of freedom. This advantage manifests itself in the degree of mesh 
refinement required to achieve a given level of accuracy. Whereas a 
FINE or EXTRA FINE mesh was required to achieve acceptable results using 
CQUA02 elements (c.f. reference 17), a COARSE or MEDIUM mesh of CQUAD4 
elements is sufficient. Such a typical mesh is shown in Figures 10 and 
11 for a T-Joint and a K-Joint respectively. The mesh was automatically 
generated using the TK00INT program described in reference (17). That 
program has recently been recoded to allow generation of the appropriate 
SHLNORM bulk data cards at each substructure grid point. 
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It should be pointed out that the CQUAD4 element does not specifi- 
cally address the problem at the intersection line of the tubular mem- 
bers. Here the question is not whether the local behavior is more 
closely approximated by thin-shell or thick-’ toil theory, but how best 
to provide a transition to the three-dimension^ 1 state of stress which 
exists and how to include the weld geometry. Again the CQUAD4 element 
has an advantage over the CQUAD2 since it has been derived from a 20- 
node hexahedron. If a mesh of these 20-node elements is designed for 
the locality of the intersection, the transitional behavior between the 
CQUAD4 and the CIHEX2 should be smooth due to tne compatibility of the 
basic interpolation functions. Unfortunately the transition between the 
two clement types will still require a rather complex set of MPC's to be 
generated and this problem has not been adequately addressed at the time 
of publication. 

The second source of potential application regards the structural 
modeling of the relatively thick shell cylinders and panels which com- 
prise the base and towers of gravity structures (see Figure 12). Sec- 
tion A-A in Figure 12 illustrates the shell connections where the pres- 
ent element could be used. The individual cells are on the order of 
20 meters in diameter and from 0,5 meters to 1,0 meters thick. The 
CQUAD2 element would be incapable of accurately modeling the structural 
behavior associated with this geometry. The CQUAD4 element provides the 
possibility of coupling with the CIHEX2 element to perform global 
stress analyses of these structures. 
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APPENDIX A 


The biquadratic interpolation functions are well-known throughout 
the literature. They are repeated here along with their derivative 
forms only for the sake of completeness. 
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Corner Nodes (i = 1, 2, 3, 4) 

- 3 (i + 44 n O O + nn n *> (44-f + nn^ “ D (A-i) 

9N i /94 = \ 4 i (i + nnp (2U., + nnp (a-2) 

aN ./an = ^ U + M n O (44^ + 2nn^) (a-3) 

Midside Nodes with | . = 0 (i = 5, 7) 

N i = \ (i - 4 2 ) a + nnp (A-4) 

9N./a^ = - 4 (1 + nnp <a-5) 

BN^/sn = | Hi (1 - 4 2 ) (A-6) 

Midside Nodes with r^- =0 (1 = 6, 8) 

N 1 = f (1 + tt t ) (1 - n 2 ) (A-7) 

9N/94 - \ Cl - n 2 ) (A-8) 

8N 1 /an = - o (1 + 44^ (A-9) 
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APPENDIX B 


At each node we are given the vector V 3 .. which is normal to the 


midsurface at that point. With reference to the basic (x, y, z) coordinate 

A A 

system, let i = (1 0 0) and j = (0 1 0). 


Choose f i; = i x which 
makes V^.. perpendicular to and the x axis. If V^. is parallel to 
i, then choose = j x t?,.. to remove the ambiguity. The third vector 


then choose V^. = j x V 3 - to remove the ambiguity. 

-l . * i j _ v^r Tvl 


of the triad is then = Vg. x V.^.. 


To compute the coordinate system transformation matrix, , normal- 
ize the components of each vector by its scalar length and form the set 
[Vii v 2i v 3i ] where 


li 


* Vl* 


li' 1 li 


etc. 
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APPENDIX C 


Equations (3-6^ define the components of the Jacobian as: 


f 


ij] =U T | 

if 

The schematic form of the inverse may be written as: 


(C-l) 


[J] 


-1 _ 


[(t x n) (n x s) (s x t)] 


det [0] 

Equations (7-9) define the local transform [0] as 
[0] = [V s \] 


(C-2) 


(C-3) 


where was computed as the normal to the surface £ = const by taking 

s x t/|s x l\ and v as well as were defined to be perpendicular to 

v . It is therefore clear that v and v. will lie in the same plane as 
-> n * .> s t r 

s and t but that v n may not in general be considered parallel to n. 

T -1 

Consider the computation of [4>] = [0] [J ] 


^ det 



[(t x n) (n x s) (s x £)] 


Now we know that (s x I) = v = |v 


5 

V 

% T 


-> 

V 


(n x s) 
(n x s) 
(n x s) 


V 

V 


(C-4) 


(s x t) 

(5 X t) 
(S X t) 


(c- 


Therefore, since the dot 
-> T 


n ' n -> | -> 

product of perpendicular vectors is zero, we have v g * (s x t) s 0. 

This completes the derivation of [<)>]. Notice that the terms and 

if’.,, are not set to zero as was done in reference (4). The only time 

that these terms would be zero is when the vector n is exactly normal to 
the surface at the point (£, q, £). This event will only occur in the 
case of flat plates. The consequence of neglecting these two terms is 
to introduce an imbalance in the moment equilibrium of the shell (c.f. 
ref. 18). 



APPENDIX D 


We choose to divide the [B 1 ] matrix into the following nodal 
partitions: 


e' * [B' B' . . . 



CD-I) 


We shall write out the expression for a typical partition [Bi] by 
rearranging appropriate terms from equation (17): 


[Bi] = 


t.. 


iT ■*: i 


t. 


Define 

g li 

- 4*^ 9N./94 + 

^2 9N i /a n 



g 2i 

= <(» 21 3N./34 + 

((>22 ®N./8r) 



g 3i 

= ^31 al V a £ + 

4>32 aN./an 



h 3i 

= ^33 N i 

hi 0 

O' ~ 




c g 21 

0 

So 

[G^] 

= m n. = 

hi hi 

0 




9 3i 0 

g li 




o g 3 , 

g 2i 



CD-3) 


CD-4) 


CD-2) 
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and 




0 

0 

0 


h 3i 


3i 


0 

0 

0 

0 

0 


(0-5) 


Fi nal ly , 


IB\) = 


t. 


[G.] re] 1 + Y & [G il + t H -jl ) [ 0 ] T V - Y (£ [G i ] + [H. ] ) [0] 


t. 
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APPENDIX E 


BULK DATA DECK 


Input Data Card CQUAD4 Quadrilateral Element Connection 


Description : Defines a homogeneous quadrilateral membrane and bending 

element (QUAD4) of the structural model. 


Format and Example : 


1 2 3 4 56789 10 


CQUAD4 

EID 

MID 

G1 

G2 

G3 

G4 

G5 

G6 

abc 

CQUAD4 

72 

13 

13 

14 

15 

16 

21 

22 

ABC 

+bc 

G7 

G8 

TH 





1 

1 


+BC 

23 

i 

24 

29.2 





i 



Field 

EID 

MID 


Content s 

Element identification number (Integer > 0) 

Identification number of a MAT1 material card (Default is 
EID) (Integer > 0) 


G1,G2,G3,G4 Grid point identification numbers of connection points 

G5 ,G6 ,G7 ,G8 (Integer > 0; G1 * G2 4 G3 t G4 * G5 * G6 * G7 ^ G8) 


TH 


Material property orientation angle in degrees (Real) 
The sketch below gives the sign convention for TH. 
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R emarks : 

1. Element identification numbers must be unique with respect to al 1 
other element identification numbers. 

2. Grid points G1 through G4 are corner nodes and must be ordered 
consecutively around the perimeter of the element in a counter 
clockwise direction. G5 through GS are midside nodes and must have 
similar ordering where: 

G5 lies between G1 and G2 
G6 lies between G2 and G3 
G7 lies between G3 and G4 
G8 lies between G4 and G1 

3. The continuation card must be present. 
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BULK DATA DECK 


Input Data Card PL0AD4 Pressure Load 


Description : Defines a uniform static pressure load applied to two- 
dimensional elements. Only QUAD4 elements may have a pressure load 
applied to them via this card. 


Format and Example : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

PL0AD4 

SID 

P 

EID 

EID 

EID 

EID 

EID 

EID 


PL0AD4 

21 

-3.6 


4 

16 


2 




Alternate Form 


PL0AD4 

SID 

P 

EID1 

"THRU" 

EID2 





PL0AD4 

1 

30.4 

16 

THRU 

48 


! — . 




Field Contents 

SID Load set identification number (Integer > 0) 

P Pressure value (Real), positive pressure value indicates 

pressure in the negative normal direction. 

EID 

EID1 Element identification number (Integer > 0; EID1 < EID2) 

EID2 


Remarks : 

1. EID must be 0 or blank for omitted entrys. 

2. Load sets must be selected in the Case Control Deck (L0AD=SID) to 
be used by NASTRAN. 

3. At least one positive EID must be present on each PL0AD4 card. 

4. If 1 he alternate form is used, all elements in the range EID1 
through EID2 must be present. 


174 





5. The "work equivalent" load vector is computed for each element 
using the relation 


a ^ ^ - T ■* 

? = P J f [N ] 1 v det [J] d| d n 

-I -1 n 


6. All elements referenced must exist. 
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BULK DATA DECK 


Input Data Card 5HLN0RM Shell Normal 


Description : Defines the direction of a normal to the shell of the 

structural model . 


Format and Example : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

SHLNORM 

ID 

CP 

XI 

X2 

X3 





SHLNORM 

2 

3 

1.0 

2.0 

3.0 






Field Contents 

ID Grid point identification number (0 < Integer < 999999) 

at which this normal is located. 


CP Identification number of coordinate system in which the 

shell normal is defined (Integer > 0 or blank ). 

X1,X2,X3 Components of the shell normal in coordinate system CP 

(Real). 


Remarks : 

1. All grid point identification numbers must be unique with respect 
to al 1 other structural, scalar, and fluid points. 

2. The meaning of XI, X2 and X3 depend on the type of coordinate 

system, CP, as follows: (see C0RD card descriptions). 


Type 

XI 

X2 

X3 

Rectangular 

X 

Y 

Z 

r ylindrical 

R 

©(degrees) 

Z 

Spherical 

R 

©(degrees) 

«t>( degrees) 
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TABLE 1, DISPLACEMENT CONVERGENCE FOR CYLINDRICAL SHELL ROOF 


GRID 

V A 

w B 

U B 

w c 


1 x 1 

-0.226 

-9.609 

-5.413 

1.801 

cm 


-0.089 

-3.783 

-2.131 

0.709 

in. 

■ 

-0.368 

-8.966 

-4.752 

1.313 

cm 

■■ 

-0.145 

-3.530 

-1.871 

0.517 

in. 

3x3 

-0.381 

-9.241 

-4.279 

1.346 

cm 


-0.150 

-3.638 

-1.921 

0.530 

in. 

4x4 

-0.381 

-9.208 

-4.854 

1.379 

cm 


-0.150 

-3.625 

-1.911 

0.543 

in. 

EXACT 

-0.384 

-9.406 

-4.986 

1.334 

cm 


-0.150 

| -3.703 

-1.963 

0.525 

in. 
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FIGURE 1. ELEMENT GEOMETRY AND LOCAL COORDINATE SYSTEMS 




FIGURE 2. CONVERGENCE STUDY FOR SIMPLY-SUPPORTED FLAT PLATE 



E = 69.4 KN/mm 2 (1.0 x 10 7 psi) 
i> = 0.3 

a = b = 254 mm (10 in.) 
t = 2.54 mm {0.1 in.) 

P = 178.4 N (40 lb.}, 6.94 KN/m 2 (1.0 psi) 



VERTICAL DEFLECTION AT CENTER 

GRID 

CENTRAL LOAD 

UNIFORM PRESSURE 


mm 

in. 

mm 

in. 

1 x 1 

- 1.286 

-5.062 x 10‘ 2 

-1.214 

-4.778 x 10' 2 

2x2 

-1.284 

-5.054 x 10' 2 

-1.142 

-4.498 x llr 2 ) 

3x3 

-1.294 

-5.095 x 10* 2 

-1.743 

-4.50 x 10' 2 

4x4 

-1.295 

-5.098 x 10' 2 

-1.139 

-4.483 x 10‘ 2 

TIMOSHENKO 

-1.287 

-5.068 x 10~ 2 

-1.128 

-4.44 x 10' 2 


ORIGINAL PAGE IS 

ot' poor quality 
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FIGURE 3. CONVERGENCE STUDY FOR PINCHED CYLINDER PROBLEM 



E « 72.9 KN/mm 2 (10.5 x 10 6 psi) 
v = 0.3125 

R = 125.8 mm (4.953 in.} 

L = 262.9 mm (10.35 in.) 
t = 2.39 mm (0.094 in.) 

P = 445.9 N (100 lb.) 


GRID 

DOF 

(TOTAL) 

RADIAL DISPLACEMENT 
AT POINT D 

mm 

in. 

1 x 1 

40 

-2.517 

-0.991 x 10' 1 

2x2 

105 

-Z753 

-1.084 x 10' 1 

3x3 

200 

-2.863 

-1.127 x 10' 1 

4x4 

325 

-2.893 

-1.139 x 10' 1 

TIMOSHENKO 


-2.761 

i 

-1.087 x 10* 1 






ORIGINAL RAGS © 
OF POOR QUALIFY 



FIGURE S. ELEMENT CONVERGENCE COMPARISON 
FOR CYLINDRICAL ROOF PROBLEM 


in 


3.5 \- 


VERTICAL 
DISPLACEMENT 
AT CENTER OF 
FREE EDGE 


3.0 


2.5 h 



cm 



TOTAL DEGREES OF FREEDOM 
FOR % MODEL 



FIGURE 6. FREE VIBRATION OF A CANTILEVER PLATE 



E = 208.3 KN/mrtr (3.0 x 10 7 psi) 
v - 0.3 

e = 7.85 t/m 3 (0.283 lb/in 3 ) 

L = 50.8 mm (2.0 in.) 
b = 25.4 mm (1.0 in.) 
x = 2,54 mm (0.1 in,} 



D/ptL fl 


EXPERIMENTAL 

NONCONFORMING TRIANGLE 



MODE 

(PLUNKETT) 

ZIEN KIEWICZ 

CQUAD4 ELEMENT 



2 x 1 

4x2 

2 x 1 

■OH 

1 

3.50 

3.39 

3.44 

3.45 

3.43 

2 

14.50 

15.30 

14,76 

14.64 

14.43 

3 

21.70 

21.16 

21.60 

22.63 

21.30 

A \ 

48.10 

49.47 

48.28 

49.79 

46.82 

5 

60.50 

67.46 

60,56 


60.55 

6 

92.30 


88.84 


90.76 

7 

92.80 


92.24 


97.17 

8 

118.70 


117.72 


123.05 

9 

125.10 i 


118.96 


130.23 

10 

154.00 






9 'Z 

s»? 

d 


















FIGURE 7. DEGENERACY RESULTING FROM NON-RECTANGULARITY 



E = 69.4 KN/mm 2 (1.0 x 10 7 psi) 
v = 0.3 

a = b = 254 mm (10 in.} 
t = 2.54 mm {0.1 in.) 

P = 6.94 KN/m 2 (1.0 psi) 


GRID 

VERTICAL DISPLACEMENT 
AT POINT D 

mm 

in. 

RECTANGULAR 

-1.143 

-4.498 x 10' 2 

10° OFFSET 

-1.145 

-4.506 x 10' 2 

20° OFFSET 

-1.118 

-4.403 x 10’ 2 

30° OFFSET 

-1.059 

-4.169 x 10' 2 

TIMOSHENKO 

-1.128 

’ -4.44 x 10" 2 
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FIGURE 8. DEGENERACY DUE TO THINNESS RATIO 


E = 69.4 KN/mm 2 fl.O x 10 7 psi) 
v = 0.3 

R = 25.4 cm (10.0 in.) 

L = 50.8 cm (20.0 in.) 

P = 445.9 N (100 lb) 


~ ~ (0.149) P ; r3 -' q ax D = - 1.0715 x 10 6 N mm 2 

2 2D (L/2 ^max 



t/R 

D 

(N mm) 

^max 

(mm) 


0.1 

« 

1.041 x 10 B 

-1.185 x 10' 2 

-1.234 x 10 s 

0.01 

1.041 x 10 5 j 

-1.035 x IQ 1 

-1.078 x 10 6 

0.001 

1.041 x 10 2 

-1.024 x 10 4 

-1.066 x 10 6 

0.0001 

1.041 x 10' 1 

-1.024 x 10 7 

L — 

-1.066 x 10 6 
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FIGURE 9. DEGENERACY L<JE TO MISALIGNMENT 
WITH DIRECTIONS OF PRINCIPAL CURVATURE 



*NOTE USE OF DEVELOPED SURFACE COORDINATES 
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FIGURE 10. SUBSTRUCTURE MESHES FOR TYPICAL T-JOINT 



ORIGINAL PAGE IS 
OF POOR QUALITY 



FIGURE 11. SUBSTRUCTURE MESHES FOR TYPICAL K-JOINT 




CHORD 


NOTE USE OF DEVELOPED SURFACE COORDINATES 



FIGURE 12. TYPICAL CONCRETE GRAVITY STRUCTURE 



')Od oW 




